clc;clear
%% Constants
G           = 6.674e-11 * ((1/1000)^3);
mass_Mars   = 6.39e23;
mass_Deimos = 1.8e15;
mass_Phobos = 10.8e15;
GM_Mars     = G*mass_Mars;
GM_Deimos   = G*mass_Deimos;
GM_Phobos   = G*mass_Phobos;
GM_sun      = 1.327124400*10^(11);
%[km^3/s^2] Gravitational Parameters

u_sm = GM_Mars / (GM_Mars+GM_sun);
u_md = GM_Deimos / (GM_Deimos+GM_Mars);
u_mp = GM_Phobos / (GM_Phobos+GM_Mars);
%[] Mu constant
%    sm --> Sun-Mars system
%    em --> Earth-Moon system%%

DU_sm = 235.34e6;
DU_mp = (9236.53+9517.58)/2;
DU_md = (23455.5+23470.9)/2;
%[km] Distance from Sun to Mars

TU_sm = sqrt(DU_sm^3/GM_sun);
TU_mp = sqrt(DU_mp^3/GM_Mars);
TU_md = sqrt(DU_md^3/GM_Mars);
%[s] Time unit

epsilon=2*10^(-8);
%[] Parameter for altering the halo orbit's state



%%
a = importdata('PhobosData.txt');
time_jd_i = str2double(char(a.textdata{1,1}));
time_utc = CalendarDate(time_jd_i)
time_utc_f = [2020,1,1,00,00,01];
time_jd_f = JulianDate(time_utc_f);
t_vec = linspace(time_jd_i,time_jd_f,3000);
integrationTime = zeros(length(t_vec),1);
for ii = 1:length(t_vec)
    
    temp = etime(CalendarDate(t_vec(ii)),time_utc);
    integrationTime(ii) = temp;
    
end

S_mp = a.data(1,:);
COE = State2Coe(S_mp',GM_Mars);
COE_circularized_phobos = [DU_mp,0,0,COE(4:6)']';
options = odeset('RelTol',1e-15,'AbsTol',1e-15);
[~,S] = ode45(@(t,S)R2bpCartesianModel(t,S,GM_Mars),integrationTime,Coe2State(COE_circularized_phobos,GM_Mars),options);

S = S';
Rpm = S(1:3,:);
Vpm = S(4:6,:);
JD_pm = t_vec;
save('PhobosEphemerisData.mat','Rpm','Vpm','JD_pm');
%%
a = importdata('DeimosData.txt');
time_jd_i = str2double(char(a.textdata{1,1}));
time_utc = CalendarDate(time_jd_i)
time_utc_f = [2020,1,1,00,00,01];
time_jd_f = JulianDate(time_utc_f);
t_vec = linspace(time_jd_i,time_jd_f,3000);
integrationTime = zeros(length(t_vec),1);
for ii = 1:length(t_vec)
    
    temp = etime(CalendarDate(t_vec(ii)),time_utc);
    integrationTime(ii) = temp;
    
end

S_mp = a.data(1,:);
COE = State2Coe(S_mp',GM_Mars);
COE_circularized_deimos = [DU_md,0,0,COE(4:6)']';
options = odeset('RelTol',1e-15,'AbsTol',1e-15);
[~,S] = ode45(@(t,S)R2bpCartesianModel(t,S,GM_Mars),integrationTime,Coe2State(COE_circularized_deimos,GM_Mars),options);

S = S';
Rdm = S(1:3,:);
Vdm = S(4:6,:);
JD_dm = t_vec;
save('DeimosEphemerisData.mat','Rdm','Vdm','JD_dm');